function [L, U, K] = periodic_tridiag_prepare(A)
% A is a tri-diagonal sparse matrix

nx = length(A);
B = A;
B(nx,1) = 0;
B(1,nx) = 0;
[L, U] = lu(B);

G = zeros(nx, 2);
H = zeros(nx, 2);
G(nx, 1) = A( nx,1 );
G(1, 2) = A(1, nx);
H(1, 1) = 1;
H(nx, 2) = 1;

K = zeros(nx, 2);
K(:,1) = U\(L\G(:,1));
K(:,2) = U\(L\G(:,2));
K = K * inv( eye(2) + transpose(H) * K );


